www.gusucode.com > matlab编程NSCT分解 图像融合 各个融合指标评价体系 分解源码程序 > matlab编程NSCT分解 图像融合 各个融合指标评价体系 分解源码程序/NSCT/histmatch.m

    function varargout = histmatch(image,matchhist,mapl)
%    HISTMATCH    直方图规定化函数
%      matchhist为规定直方图，规定图像的灰度级N<=M
%      没有该参数输入，默认为执行直翻图均衡化
%      mapl只有有matchhist参数时才具有意义，为映射规则
%      可以取值为SML和GML，默认为SML
%      无输出默认输出原始图像、原始直方图、规定图像、规定图像直方图
%    $Author: lskyp    $Date: 2009.08.18
error(nargchk(1,3,nargin));
% 没有规定直方图输入，即输入参数只有一个，执行直方图均衡化
if nargin == 1
    if nargout == 0
        histeq_my(image);
    else
        [imageeq,grayx,histeq,hist] = histeq_my(image);
        varargout = {imageeq,grayx,histeq,hist};
    end
else
    if nargin == 2
        mapl = 'GML';
    end
    if ~ (strcmp(mapl,'GML') || strcmp(mapl,'SML'))
        error('错误的映射规则')
    end
    if ndims(matchhist) > 2 || (ndims(matchhist) == 2 && (size(matchhist,1) ...
            >1 && size(matchhist,2) > 1))
        error('错误的规定直方图，规定直方图需为向量')
    end
    matchhist = matchhist(:);
    matchhist = matchhist';
    N = length(matchhist);
    M = double(intmax(class(image))) + 1;
    if N > M
        error('错误的规定直方图，规定直方图灰度级需<=原始灰度级')
    elseif N < M % 可以视为>N的灰度部分置零
        matchhist = [matchhist zeros(1,M - N)];
    end
    [hist,grayx] = imhist_my(image);
    % 归一化
    histn = hist/numel(image);
    % 原始直方图的累积直方图
    histc = cumsum(histn);
    % 规定直方图的累积直方图
    matchhistc = cumsum(matchhist);
    % 映射的实现，Ref. 章毓晋. 图像工程（上册）——图像处理
    tk = zeros(1,M); % 映射关系hist-->tk
    if strcmp(mapl,'SML')
        tk(1) = 1; % 2009.08.19增加，为下一步修改准备
        linit = 2;
        kinit = 0;
        for k = 1:M
            if histc(k) == 0
                kinit = kinit + 1;
            else
                break;
            end
        end
        
        % 标记为SML为零的位置
        indis0 = zeros(1,N);
        
        for k = kinit:M
            mappingL = zeros(size(matchhist));
            % 2009.08.19为节省运算量而修改
            for l = linit:N % 2009.08.19修改，规定直方图从上一次最小处计算，节省运算时间
                if matchhistc(l) == matchhistc(l - 1);
                    continue;
                end
                mappingL(l) = abs(histc(k) - matchhistc(l));
                if mappingL(l) == 0
                    indis0(linit) = 1;
                end
            end
            mappingL(N) = abs(histc(k) - matchhistc(N));
            indtemp1 = find(mappingL == 0);
            indtemp2 = find(indis0 == 1);
            mappingL(setdiff(indtemp1,indtemp2)) = inf;
            [minmap,tk(k)] = min(mappingL(linit:N));
            linit = tk(k) + linit - 1;
            tk(k) = linit;
            for k = 1:kinit - 1
                tk(k) = tk(kinit);
            end
        end
    else % 组映射规则GML
        % 2009.08.19 Ref. 章毓晋. 《图像工程》（上册） P96，图4.4.8 GML映射计算图
        graymarker = 1;
        for k = 1:N
            if matchhistc(k) == 0;
                graymarker = graymarker + 1;
            else
                break;
            end
        end
        linit = graymarker;
        tk(1) = linit;
        tkcnt = 2;
        
        % 标记为GML为零的位置
        indis0 = zeros(1,N);
        
        for l = linit:N
            mappingL = zeros(1,M);
            if matchhistc(l) == matchhistc(l - 1);
                continue;
            end
            for k = 1:M - 1
                if histc(k) == histc(k + 1) % 跳过原始累积直方图有相同值的元素
                    % 如果不跳过，那么再寻找最小值时，会有一系列最小值，而MATLAB
                    % 只找到第一个最小值的位置，但是需要的应该是最后一个，跳过后，
                    % mappingL矩阵中会出现0值，那么寻找最小值会误算，因此，在下面
                    % 的语句中有一句将这些零值赋为最大inf
                    continue;
                end
                mappingL(k) = abs(matchhistc(l) - histc(k)); % GML规则
                if mappingL(k) == 0
                    indis0(k) = 0;
                end
            end
            mappingL(M) = abs(matchhistc(l) - histc(M));
            indtemp1 = find(mappingL == 0);
            indtemp2 = find(indis0 == 1);
            mappingL(setdiff(indtemp1,indtemp2)) = inf;
            [minmap,minind] = min(mappingL);
            for pp = tkcnt:minind
                tk(pp) = l; % 建立灰度映射矩阵
            end
            tkcnt = minind + 1;
        end
    end
    if N < M
        for k = N:M
            tk(k) = N;
        end
    end
    histm = zeros(1,N);
    for k = 1:M
        histm(tk(k)) = histm(tk(k)) + hist(k);
    end
    tk = tk - 1; % 方便灰度值的计算
    % 图像的规定化
    imagematch = zeros(size(image));
    for p = 1:size(image,1)
        for q = 1:size(image,2)
            ind = image(p,q) + 1;
            imagematch(p,q) = tk(ind);
        end
    end
    eval(['imagematch = ' class(image) '(imagematch);'])
    graymatch = 0:length(histm) - 1;
    % 输出
    if nargout == 0
        subplot(221)
        imshow(image)
        title('原始图像')
        subplot(222)
        imshow(imagematch)
        title('直方图规定化后图像')
        subplot(223)
        stem(grayx,hist,'Marker','none')
        xlim([0 max(grayx)]);
        xlabel('灰度值')
        ylabel('统计个数');
        title('原始直方图')
        subplot(224)
        stem(graymatch,histm,'Marker','none')
        xlim([0 max(graymatch)]);
        xlabel('灰度值')
        ylabel('统计个数');
        title(['规定化后直方图' mapl])
    else
        varargout = {imagematch,graymatch,histm,grayx,hist};
    end
end